Hysteresis in one-dimensional reaction-diffusion systems 
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We introduce a simple nonequilibrium model for a driven diffusive system with nonconservative 
reaction kinetics which exhibits ergodicity breaking and hysteresis in one dimension. These phe- 
nomena can be understood through a description of the dominant stochastic many-body dynamics 
in terms of an equilibrium single-particle problem, viz. the random motion of a shock in an effective 
potential. This picture also leads to the exact phase diagram of the system and suggests a new 
generic mechanism for "freezing by heating" . 
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The closely related questions of phase coexistence, er- 
godicity breaking and hysteresis in noisy one-dimensional 
systems with short range interactions and finite local 
state space (such as in spin systems or vertex models) are 
intriguing. In thermal equilibrium these phenomena can- 
not occur as there is no local mechanism that could limit 
the growth of islands of a minority phase inside a ma- 
jority phase. Far from equilibrium one has found phase 
separation and spontaneous symmetry breaking in driven 
diffusive systems provided that either a bulk conservation 
law, viz. particle number conservation 0, H H LH Q , or 
vanishing local transition rates [a, L| constrain the local 
dynamics. As already noted in Ref. 0] the only known 
exception to this rule, the error-correcting model by Gacs 
[U , is rather complicated and still not widely understood, 
see also 0- 

Recently it has been demonstrated that phase coex- 
istence occurs in a one-dimensional driven diffusive sys- 
tem even in the presence of Langmuir kinetics A ^ 
which break the bulk conservation law [loj . This mecha- 
nism is inspired by the process of motor proteins moving 
along actin filaments. Earlier this model was introduced 
as a toy model reproducing stylized facts in limit order 
markets jjjj . The formation of a localized shock in this 
system which separates a domain of low particle density 
from a domain of high density has been studied subse- 
quently However, the two different domains do 
not represent two possible global steady states. The pro- 
cess is ergodic even in the thermodynamic limit and no 
hysteresis is possible. 

It is the purpose of this letter to present a simple 
nonequilibrium system with local non-conservative dy- 
namics and finite local state space which exhibits er- 
godicity breaking and hysteresis in the thermodynamic 
limit, in the usual sense that in finite volume the so- 
journ time in two metastable steady states increases ex- 
ponentially with system size. To be specific we investi- 
gate the totally asymmetric exclusion process (TASEP) 
augmented by nonconservative reaction kinetics. The 
TASEP is a stochastic model of diffusing particles on 
a one-dimensional lattice with a hopping bias in one di- 
rection E3- Each site from 1 to L is either empty or 



occupied by one particle. In the bulk particles ('A') hop 
stochastically from site i to i + 1 with unit rate, provided 
that the target site is empty. The boundaries act as par- 
ticle reservoirs with densities p_ on the left resp. p+ on 
the right: On site 1 particles are created with rate p_, 
provided the site is empty, which corresponds to a par- 
ticle hopping from the left reservoir onto the first site. 
Particles on site L are annihilated with rate 1 — p+, cor- 
responding to a particle hopping from the last site into 
the right reservoir. 

In our model particles also undergo the following re- 
action process: On a vacant site enclosed by two par- 
ticles a particle may be attached with rate u> a , and a 
particle enclosed by two other particles may be detached 
with rate uJd- This process can be symbolically written as 
AszA f± AAA and may be interpreted as activated Lang- 
muir kinetics. Without the TASEP dynamics the station- 
ary density of this process is either K = Lo a /(Lo a + u>d) 
or zero, with no correlations |l5j| . As in previous work 
we consider the physically interesting case when L — > oo 
and these rates are proportional to l/L El El El El • 
Hence we define renormalized rates 

u> a = Q. a /L, uj d = fl d /L (1) 

where Q a and fid are kept constant while L — > oo. For 
other choices of the attachment/detachment (AD) rates 
the dynamics is either governed by the TASEP (u> a ,d < 
0(1/ L)) or by the AD process (cj m > 0(1/ L)). 

We find a stationary phase diagram of the model with 
five distinct phases (Fig. . The stationary density pro- 
file pi is not constant as a function of lattice site i. Yet 
some of the phases are comparable to those of the usual 
TASEP with open boundaries EE Ell : m the high den- 
sity phases (HD1/2) one finds pi > 1/2 while in the low 
density phase (LD) p. L < 1/2. In HD1 the bulk density 
profile is dependent on p + , while it is independent of 
both boundaries in HD2 as in the maximal current phase 
of the TASEP. On the other hand two additional phases 
exist: (i) A coexistence phase which is characterized by 
a stable shock position, i.e., a jump in the density pro- 
file which is localized at a certain position in the bulk of 
the system. The shock connects a low density domain to 
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FIG. 1: Phase diagram for £l a = 0.7 and Qd = 0.1 with two 
high density phases (HDI, HD2), a low density phase (LD), 
a coexistence phase and the nonergodic phase. 



0.8 




P_ 

FIG. 2: Hysteresis plot for L = 2000, Q a = 0.7, D. d = 0.1, 
p+ = 0.45. p- was changed by 10 -4 in every 5000 (solid 
line), 1500 (dashed line) and 500 (dotted line) MC steps. The 
hysteresis loop gets wider as the speed of changing p_ is in- 
creased. 



its left with a high density domain to its rig ht as known 
from related models studied previously |l0|, |l2| . Notice 
that in the usual TASEP there is a coexistence line in the 
phase diagram with a nonlocalized shock. In a different 
parameter regime we find a novel phase with an unstable 
shock position in the bulk. In this phase both the LD and 
HD states are stable (if L — > oo) which implies that er- 
godicity is broken in the thermodynamic limit. Although 
for finite systems a transition between the two states is 
possible, the mean life time of each steady state is ex- 
ponentially large in the system size L (see below). We 
note that this is not a spontaneous symmetry breaking 
since there is no symmetry relating the two metastable 
states. This phase has no analog in the TASEP with 
open boundaries. 

Hysteresis in this nonequilibrium setting was observed 
by measuring the space-averaged density p along the 
curve of constant p + = 0.45 while changing p_ in such a 
way that the system starting from the LD phase passed 
through the nonergodic phase and ended up in the HD2 
phase. Then the process of changing p_ was reversed. A 
relevant parameter in hysteresis phenomena is the speed 
of sweeping: in our simulations p_ was changed by 10~ 4 
in every k MC steps (k = 500,1500,5000). A time av- 
erage was not taken, p was measured in every k steps. 
On Fig. [21 one can see the resulting hysteresis loops. We 
found that the hysteresis loop inflates with increasing 
speed which is reminiscent of hysteresis in usual mag- 
netic systems. 

To rationalize these observations we first consider the 
hydrodynamic limit on the Euler scale, i.e., we take L — > 
oo while the lattice constant is scaled by a = 1/L and the 
time by t = t\ ai f lcc / L . Thus the spatial coordinate x = 
i/L becomes continuous. Following the line of arguments 
of Ref. |l2j the hydrodynamic equation for the density 



takes the form 



d d 



(2) 



with the exact current j(p) = p(l — p) of the TASEP and 
the cubic source term 



S( P ) = n aP 2 {i - P ) - n dP 3 . 



(3) 



resulting from the activated Langmuir kinetics. For the 
stationary state dtp(x,t) = holds and using d x j = 
dj I dp ■ dp/dx we obtain 



v c (p) 



dx 



S(p), 



(4) 



with the collective velocity v c = dj /dp. This nonlinear 
differential equation can be integrated analytically and 
yields the flow field 



Xip) = ^P + ^^ lli 



1 

K 



(5) 



with an integration constant c. 

As the differential equation is of first order and the 
boundary condition fixes the density at two positions, 
following a line of the flow field does not represent a so- 
lution of the boundary problem in general. In the original 
lattice model this inconsistency is resolved by the appear- 
ance of a shock and/or boundary layers as described in 
|l2l |. Apart from the discontinuities the stationary den- 
sity profile follows the flow field of eq. 

In order to understand quantitatively the selection 
of the stationary shock position (which determines the 
phase diagram) and also to explain the phenomenon of 
hysteresis from a microscopic viewpoint we describe the 
dominant dynamical mode of the particle system in terms 
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of the random motion of the shock. To this end we gener- 
alize the approach of |l8| and introduce space-dependent 
hopping rates 



^x — *x-\-a 



3r{x) 



w 



x-\-a — *x 



Pr{x) - p L (x) 
Pr{x)- p L {x)' 



(6) 



for jumps of the shock over a lattice constant a. Here 
the indices L and R denote the solutions (lines of the 
flow field 10) on the left resp. right of the shock. Similar 
hopping rates are used in \\?\ . The space-dependent hop- 
ping rates furnish us with the picture of a random walker 
in an effective energy landscape E{x) inside a finite box. 
The energy landscape is generated by the interplay of the 
particle current with the reaction kinetics. In this way we 
relate the original nonequilibrium many-particle system 
to an equilibrium single-particle model. Let p{x) be the 
equilibrium probability of the shock being at position x. 
Then due to detailed balance 



l^x— *x+a 

_ p(x + a) 

W x+a -> x P(x) 



= exp(-E(x + a)+E(x)). (7) 



which defines the energy landscape. 

The potential E{x) is monotonically increasing (de- 
creasing) function for the HD (LD) phase (Fig. 03). This 
implies that although there are fluctuations the shock is 
always driven to the left (right) boundary. In the co- 
existence phase there is a global minimum in the bulk 
resulting in a stable shock position (Fig. |3J) at a macro- 
scopic distance from the boundaries. Here the dynamics 
can be well approximated by a random walker in a har- 
monic potential which gives a Gaussian distribution for 
the shock position. Hence the width of the shock distri- 
bution is proportional to \fL 0] which was also found 
in [H El for the TASEP with Langmuir kinetics. 

The nonergodic phase is characterized by a global en- 
ergy maximum in the bulk (Fig.[3J), leading to an unstable 
bulk fixed point of the shock. The two minima at the left 
and right boundary correspond to the two stable station- 
ary states. Starting with an initial condition close to one 
of the minima, the random walker will drift most likely 
into this local minimum and stay in its vicinity for a time 
of the order of the mean first passage time f before it tra- 
verses to the other minimum. This leads to hysteresis. 
Using a formula for the mean first passage time derived 
by Murthy and Kehr one expects that f grows expo- 
nentially with the system size L. Moreover, one expects 
the transition from one minimum to the other to be a 
random Poisson process with an average waiting time f . 

This simple one-particle picture is well borne out by 
MC simulations. For judiciously chosen parameters it is 
possible to perform simulations up to times much larger 
than f . Using multispin coding [20j for the MC algorithm 
rather good statistics become available for the waiting 
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FIG. 3: Examples for the energy landscape in four phases. 
Note that in the HD and LD phases E(x) can be either convex 
or concave 
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FIG. 4: Snapshot of the time evolution of the scaled position 
of the second class particle for L = 1000, p_ = 0.2705, p+ = 
0.63, Q a = 0.5, Qd — 0.1. A position of the second class par- 
ticle near the left boundary (x m 0) corresponds to the high 
density state, while a position near the right boundary (x ~ 1) 
corresponds to the low density state. 



time t (the time the system spends in one of the sta- 
tionary states before switching to the other). For tracing 
the position of the shock we use the second class particle 
technique [ilj . We measured the position of the second 
class particle as a function of time: a typical realization 
is shown in Fig. 0] 

As shown in Fig.[5]the numerically determined cumula- 
tive distribution function $(f) = P(t < t) of the waiting 
time r is hardly distinguishable from the expected expo- 
nential distribution 



$(t) = 1 - exp(-t/f). 



(8) 



With this picture of a moving shock in mind and us- 
ing the expression (0 it is also possible to derive the 
exact phase transition lines defining the phase diag ram 
presented above. Adapting the arguments of [12( the 
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FIG. 5: Numerically determined cumulative distribution of 
the transition times from the upper state to the lower (solid 
line) compared to the exponential distribution (dashed line) 
with parameters as in 2] Similar results are found for the 
transition in the other direction, but with a different f ITBI . 



absence of this noise, i.e., in the usual TASEP, the shock 
performs an unbiased random walk and hence is unlo- 
calized, whereas suitably chosen reaction kinetics may 
create a variety of effective potentials which localize the 
shock. An increase in noise strength is usually associ- 
ated with heating up a system whereas localization re- 
duces the amount of disorder, corresponding to cooling. 
Thus we have identified a novel generic mechanism for 
the phenomenon of freezing by heating. The description 
of the nonequilibrium many-body dynamics in terms of a 
collective single-particle mode moving under equilibrium 
conditions yields the exact stationary phase diagram as 
well as the numerically verified flipping process between 
the metastable states of the finite system. Details of the 
flipping dynamics will be presented elsewhere .15J. 

A. Rakos thanks the Deutsche Forschungsgemeinschaft 
(DFG) for financial support. 



analysis of the stability of boundary layers yields a high 
density phase for p_ > 1/2. Because of a boundary 
layer the bulk solution of the density profile is inde- 
pendent of p+ if p+ < 1/2. Thus in this region the 
phase diagram is only ruled by For p_ < 1/2 the 
two lines in the phase diagram bounding the coexistence 
phase and nonergodic phase resp. are determined by the 
stationary shock position. Crossing the line separating 
LD/coexistence phase and nonergodic/HD phase from 
left to right results in a change of the sign of d x E(l) 
from — — * +. Crossing the other line separating co- 
existence/HD phase and LD/nonergodic phase from left 
to right results in a change of the sign of d x E(0) from 

► +. The sign of the slope of the energy profile, i.e., 

the stability of the shock position can be analysed by 
considering the average shock velocity 



3r(x) ~ 3l{x) 
Pr{x) - Pl(x)' 



(9) 



A shock position at the boundary is stable when it is 
driven toward the boundary, i.e., v s (0) < at the left, 
v s (l) > at the right boundary. Thus the lines separat- 
ing the phases are calculated by comparing the values of 
Pl(x) and Pr(x) at the positions x = 0, 1. 

To conclude we have demonstrated the existence of 
hysteresis and broken ergodicity (in the thermodynamic 
limit) in a driven diffusive system without bulk conser- 
vation law. We stress that the two different stationary 
distributions are not ordered states in which the acti- 
vated Langmuir reaction kinetics would be dynamically 
suppressed. Surprisingly, adding noise which is on av- 
erage spatially homogeneous (the nonconservative reac- 
tion process) to a conservative spatially homogeneous 
nonequilibrium system with a nonvanishing particle cur- 
rent leads to a space- dependent effective potential which 
determines the stationary position of the shock. In the 
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